
/*  Triangulation of a sphere using spherical coordinates 
 *
 *  Spherical coordinates of a point P is (r, theta, phi)
 *     where      r = distance of P from the origin O
 *            theta = angle ZOP 
 *              phi = angle XOQ,   Q is the foot of the perpendicular
 *                                 from P on the XY plane.
 *     Actually, theta = PI/2 - latitude of P and
 *                 phi = longitude of P.
 *
 *  The correspondance with x, y, z coordinates:
 *                x = r sin(theta) cos(phi)
 *                y = r sin(theta) sin(phi)
 *                z = r cos(theta) 
 */

#include <stdio.h>
#include <math.h>

main()
{
   double          xcenter, ycenter, zcenter, radius;
   int             num_subdivisions;
   int             numpts, numtri, numedg;
   double           x, y, z;
   int             i, j, k, l, count;
   int             num_subdiv_long;
   int             s1, s2, s10, s20;
   double          theta, phi, del_theta, del_phi;
   double          PI = 3.14159265359;
   FILE            *out_file;
   char            file_name[12];


   /* Get center and radius  */

      printf("Give center & radius: ");
      scanf("%lf %lf %lf %lf", &xcenter, &ycenter, &zcenter, &radius);

   /* Number of subdivisions */

      printf("Give number of subdivisions: ");
      scanf("%d", &num_subdivisions);

   /* Output file name */

      printf("Give output file name: ");
      scanf("%s", file_name);
      out_file = fopen(file_name, "w");
      if (out_file == NULL)
	 {
	    printf("Can't open ouput file\n");
	    exit(1);
         }

   /* Compute total number of vertices and triangles  */

      numpts = 2 + 4 * num_subdivisions * num_subdivisions;
      numtri = 8 * num_subdivisions * num_subdivisions;
      numedg = 12 * num_subdivisions * num_subdivisions;
      fprintf(out_file,"%d  %d  %d \n", numpts, numtri, numedg);

   /* Compute and write the coordinates of the points */

      count = 0;

      /* First the point corresponding to (r,0,0) */
      x = xcenter; y = ycenter; z = zcenter + radius;
      fprintf(out_file, "%lf  %lf  %lf \n", x, y, z);
      count = count + 1;

      /* Compute points in each latitude level */

      del_theta = PI / (2.0 * num_subdivisions);

      for (i = 1; i <= 2 * num_subdivisions - 1; i++)
         {
            theta = del_theta * (double) i;

            if (i <= num_subdivisions) num_subdiv_long = i;
	    else  num_subdiv_long = 2 * num_subdivisions - i;
            del_phi = PI / (2.0 * num_subdiv_long);

            for (j = 0; j <= 4 * num_subdiv_long - 1; j++)
	       {
		  phi = del_phi * (double) j;

		  x = xcenter + radius * sin(theta) * cos(phi);
		  y = ycenter + radius * sin(theta) * sin(phi);
		  z = zcenter + radius * cos(theta);
                  fprintf(out_file, "%lf  %lf  %lf \n", x, y, z);
                  count = count + 1;
               }
         }
      
      /* Last point corresponding to (r,PI,0) */
      x = xcenter; y = ycenter; z = zcenter - radius;
      fprintf(out_file, "%lf  %lf  %lf \n", x, y, z);
      count = count + 1;

      /* Check the number of points  */
      if (count != numpts)
	 {
	    printf("Count is not equal to number of points\n");
	    printf("count = %d numpts = %d\n", count, numpts);
	    exit(1);
         }

   /* List the triangles  */

      /* First row  */
      fprintf(out_file, "3   1   2   3\n");
      fprintf(out_file, "3   1   3   4\n");
      fprintf(out_file, "3   1   4   5\n");
      fprintf(out_file, "3   1   5   2\n");
      count = 4;

      /* Above equator */
      s1 = 2;  s2 = 6; 
      for (i = 2; i <= num_subdivisions; i++)
         {
	    s10 = s1;  s20 = s2;
	    for (l = 1; l <= 4; l++)
	       {
                  fprintf(out_file, "3     %d  %d  %d \n", s1, s2, s2+1);
	          count = count + 1;
	          s2 = s2 + 1;

	          for (k = 2; k <= i; k++)
		     if ((l < 4) || (k < i))
		        {
                           fprintf(out_file, "3     %d  %d  %d \n", s1, s2, s1+1);
                           fprintf(out_file, "3     %d  %d  %d \n", s1+1, s2, s2+1);
	                   count = count + 2;
	                   s1 = s1 + 1; s2 = s2 + 1;
                        }			
		     else
		        {
                           fprintf(out_file, "3     %d  %d  %d \n", s1, s2, s10);
                           fprintf(out_file, "3     %d  %d  %d \n", s10, s2, s20);
	                   count = count + 2;
	                   s1 = s1 + 1; s2 = s2 + 1;
                        }
                }
           }

      /* Below equator */
      for (i = num_subdivisions; i >= 2; i--)
	 {
	    s10 = s1;  s20 = s2;
	    for (l = 1; l <= 4; l++)
	       {
                  fprintf(out_file, "3     %d  %d  %d \n", s1, s2, s1+1);
	          count = count + 1;
	          s1 = s1 + 1;
   
	          for (k = 2; k <= i; k++)
		     if ((l < 4) || (k < i))
		        {
                           fprintf(out_file, "3     %d  %d  %d \n", s1, s2, s2+1);
                           fprintf(out_file, "3     %d  %d  %d \n", s1, s2+1, s1+1);
	                   count = count + 2;
	                   s1 = s1 + 1; s2 = s2 + 1;
                        }			
		     else
		        {
                           fprintf(out_file, "3     %d  %d  %d \n", s1, s2, s20);
                           fprintf(out_file, "3     %d  %d  %d \n", s1, s20, s10);
	                   count = count + 2;
	                   s1 = s1 + 1; s2 = s2 + 1;
                        }
                }
           }

      /* Last row  */
      fprintf(out_file, "3     %d  %d  %d \n", s1, s2, s1+1);
      fprintf(out_file, "3     %d  %d  %d \n", s1+1, s2, s1+2);
      fprintf(out_file, "3     %d  %d  %d \n", s1+2, s2, s1+3);
      fprintf(out_file, "3     %d  %d  %d \n", s1+3, s2, s1);
      count = count + 4;

      /* Check the number of triangles  */
      if (count != numtri)
	 {
	    printf("Count is not equal to number of triangles\n");
	    printf("count = %d numtri = %d\n", count, numtri);
	    exit(1);
         }
}

